﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace CloudChartLib
{
    /// <summary>
    /// 插值
    /// </summary>
    public class Interpolation
    {
        /// <summary>
        /// 二维数组插值
        /// </summary>
        /// <param name="source">源二维数组</param>
        /// <param name="targetRow">目标行数</param>
        /// <param name="targetCol">目标列数</param>
        /// <param name="method">插值方法</param>
        /// <returns>目标二维数组</returns>
        public static double[,] Interpolate(double[,] source, int targetRow, int targetCol, InterpolateMethod method)
        {
            int sourceRow = source.GetLength(0);
            int sourceCol = source.GetLength(1);

            double[,] target = new double[targetRow, targetCol];

            if (method == InterpolateMethod.Nearest)
            {
                double ratioRow = (double)sourceRow / targetRow;
                double ratioCol = (double)sourceCol / targetCol;

                for (int i = 0; i < targetRow; i++)
                {
                    for (int j = 0; j < targetCol; j++)
                    {
                        target[i, j] = source[(int)(i * ratioRow), (int)(j * ratioCol)];
                    }
                }
            }
            else if (method == InterpolateMethod.Bilinear)
            {
                double ratioRow = (double)(sourceRow - 1) / (targetRow - 1);
                double ratioCol = (double)(sourceCol - 1) / (targetCol - 1);

                for (int i = 0; i < targetRow; i++)
                {
                    double sR = i * ratioRow;
                    int iR = (int)sR;
                    double dR = sR - iR;

                    for (int j = 0; j < targetCol; j++)
                    {
                        double sC = j * ratioCol;
                        int iC = (int)sC;
                        double dC = sC - iC;

                        double value = (1 - dR) * (1 - dC) * source[iR, iC];
                        if (dC != 0)
                        {
                            value += (1 - dR) * dC * source[iR, iC + 1];
                        }
                        if (dR != 0)
                        {
                            value += dR * (1 - dC) * source[iR + 1, iC];
                        }
                        if (dC != 0 && dR != 0)
                        {
                            value += dR * dC * source[iR + 1, iC + 1];
                        }

                        target[i, j] = value;
                    }
                }
            }
            else if (method == InterpolateMethod.BiCubic)
            {
                double ratioRow = (double)sourceRow / targetRow;
                double ratioCol = (double)sourceCol / targetCol;

                int rowMax = sourceRow - 1;
                int colMax = sourceCol - 1;

                double sourceMin = double.MaxValue;
                double sourceMax = double.MinValue;
                for (int i = 0; i < sourceRow; i++)
                {
                    for (int j = 0; j < sourceCol; j++)
                    {
                        sourceMin = Math.Min(sourceMin, source[i, j]);
                        sourceMax = Math.Max(sourceMax, source[i, j]);
                    }
                }

                double targetMin = double.MaxValue;
                double targetMax = double.MinValue;

                for (int i = 0; i < targetRow; i++)
                {
                    double sR = i * ratioRow - 0.5;
                    int iR = (int)sR;
                    double dR = sR - iR;

                    for (int j = 0; j < targetCol; j++)
                    {
                        double sC = j * ratioCol - 0.5;
                        int iC = (int)sC;
                        double dC = sC - iC;

                        double value = 0;

                        for (int m = -1; m < 3; m++)
                        {
                            double k1 = BiCubicInterpolator(dR - m);

                            int iR2 = iR + m;
                            if (iR2 < 0)
                                iR2 = 0;
                            if (iR2 > rowMax)
                                iR2 = rowMax;

                            for (int n = -1; n < 3; n++)
                            {
                                double k2 = k1 * BiCubicInterpolator(n - dC);

                                int iC2 = iC + n;
                                if (iC2 < 0)
                                    iC2 = 0;
                                if (iC2 > colMax)
                                    iC2 = colMax;

                                value += k2 * source[iR2, iC2];
                            }
                        }
                        target[i, j] = value;

                        targetMin = Math.Min(targetMin, value);
                        targetMax = Math.Max(targetMax, value);
                    }
                }

                double ratio = (sourceMax - sourceMin) / (targetMax - targetMin);
                for (int i = 0; i < targetRow; i++)
                {
                    for (int j = 0; j < targetCol; j++)
                    {
                        target[i, j] = ratio * (target[i, j] - targetMin) + sourceMin;
                    }
                }
            }

            return target;
        }

        /// <summary>  
        /// 双三次插值器
        /// </summary>
        private static double BiCubicInterpolator(double x)
        {
            if (x < 0)
            {
                x = -x;
            }

            double biCoef = 0;

            if (x <= 1)
            {
                biCoef = (1.5 * x - 2.5) * x * x + 1;
            }
            else if (x < 2)
            {
                biCoef = ((-0.5 * x + 2.5) * x - 4) * x + 2;
            }

            return biCoef;
        }

        /// <summary>
        /// 三次样条插值
        /// </summary>
        /// <param name="x">离散点x</param>
        /// <param name="y">离散点y</param>
        /// <param name="xs">插值点x</param>
        /// <returns>插值点y</returns>
        public static double[] Spline(double[] x, double[] y, double[] xs)
        {
            int plength = x.Length;
            double[] h = new double[plength];
            double[] f = new double[plength];
            double[] l = new double[plength];
            double[] v = new double[plength];
            double[] g = new double[plength];

            for (int i = 0; i < plength - 1; i++)
            {
                h[i] = x[i + 1] - x[i];
                f[i] = (y[i + 1] - y[i]) / h[i];
            }

            for (int i = 1; i < plength - 1; i++)
            {
                l[i] = h[i] / (h[i - 1] + h[i]);
                v[i] = h[i - 1] / (h[i - 1] + h[i]);
                g[i] = 3 * (l[i] * f[i - 1] + v[i] * f[i]);
            }

            double[] b = new double[plength];
            double[] tem = new double[plength];
            double[] m = new double[plength];
            double f0 = (y[0] - y[1]) / (x[0] - x[1]);
            double fn = (y[plength - 1] - y[plength - 2]) / (x[plength - 1] - x[plength - 2]);

            b[1] = v[1] / 2;
            for (int i = 2; i < plength - 2; i++)
            {
                b[i] = v[i] / (2 - b[i - 1] * l[i]);
            }
            tem[1] = g[1] / 2;
            for (int i = 2; i < plength - 1; i++)
            {
                tem[i] = (g[i] - l[i] * tem[i - 1]) / (2 - l[i] * b[i - 1]);
            }
            m[plength - 2] = tem[plength - 2];
            for (int i = plength - 3; i > 0; i--)
            {
                m[i] = tem[i] - b[i] * m[i + 1];
            }
            m[0] = 3 * f[0] / 2.0;
            m[plength - 1] = fn;
            int xlength = xs.Length;
            double[] insertRes = new double[xlength];
            for (int i = 0; i < xlength; i++)
            {
                int j = 0;
                for (j = 0; j < plength; j++)
                {
                    if (xs[i] < x[j])
                        break;
                }
                j = j - 1;
                if (j == -1 || j == x.Length - 1)
                {
                    if (j == -1)
                        throw new Exception("插值下边界超出");
                    if (j == x.Length - 1 && xs[i] == x[j])
                        insertRes[i] = y[j];
                    else
                        throw new Exception("插值下边界超出");
                }
                else
                {
                    double p1;
                    p1 = (xs[i] - x[j + 1]) / (x[j] - x[j + 1]);
                    p1 = p1 * p1;
                    double p2; p2 = (xs[i] - x[j]) / (x[j + 1] - x[j]);
                    p2 = p2 * p2;
                    double p3; p3 = p1 * (1 + 2 * (xs[i] - x[j]) / (x[j + 1] - x[j])) * y[j] + p2 * (1 + 2 * (xs[i] - x[j + 1]) / (x[j] - x[j + 1])) * y[j + 1];

                    double p4; p4 = p1 * (xs[i] - x[j]) * m[j] + p2 * (xs[i] - x[j + 1]) * m[j + 1];
                    p4 = p4 + p3;
                    insertRes[i] = p4;
                }
            }
            return insertRes;
        }
    }
}
